# ================================================#
#             Replication files for:              #
#   "The Legacy of Multi-Dimensional Conflict"    #
#         Joan Barceló and Leonid Peisakhin       #
#          Comparative Political Studies          #
#          Last updated: August 29, 2025          #
# ================================================#

# ------------------
# AGGREGATE ANALYSIS: 
# Replication script prepared to reproduce Tables 1 and 2, and Figure 2 in the main text in the manuscript.
# Replication code prepared to reproduce Appendices G, H.1, I, J, K, L.1, L.3, O.1, O.3, P.1, Q.2, 
# ------------------

# Replication files prepared using: "R version 4.2.3 (2023-03-15)"#

# ------------------
# Load packages
# ------------------
library(dplyr)
library(purrr)
library(tidyr)
library(stringr)
library(broom)
library(rlang)
library(kableExtra) 
library(ggplot2)
library(tibble)
library(knitr)
library(ivreg)

# ------------------
# Load data
# ------------------

PATH <- ''

aggregate_data <- 
  readstata13::read.dta13(paste0(PATH, "aggregate_data.dta"))

# ---------------------------------- #
# ---------------------------------- #
# ---------------------------------- #
#              MAIN TEXT             #
# ---------------------------------- #
# ---------------------------------- #
# ---------------------------------- #

# ------------------#
# ------------------#
#      Table 1      #
# ------------------#
# ------------------#

stat_row <- function(var_expr, label, data = aggregate_data) {
  v <- eval(parse(text = var_expr), envir = data)
  if (!is.numeric(v)) v <- as.numeric(v)
  
  g <- data$group
  west <- v[g == "West"]
  east <- v[g == "East"]
  
  west_m <- mean(west, na.rm = TRUE)
  west_sd <- sd(west, na.rm = TRUE)
  east_m <- mean(east, na.rm = TRUE)
  east_sd <- sd(east, na.rm = TRUE)
  
  tval <- tryCatch(
    t.test(west, east, alternative = "two.sided", var.equal = FALSE)$statistic,
    error = function(e) NA_real_
  )
  
  tibble::tibble(
    label = label,
    west_mean = west_m,
    west_sd   = west_sd,
    east_mean = east_m,
    east_sd   = east_sd,
    t_stat    = abs(unname(tval))
  )
}

rows <- tribble(
  ~label, ~expr,
  
  # 1931 Spanish national elections
  "% Vote for Catalan nationalist parties, 1931",                    "catxleft_1931 + catxright_1931",
  "% Vote for Left parties, 1931",                       "spnxleft_1931 + catxleft_1931",
  "% Vote for the Spanish Left, 1931",                   "spnxleft_1931",
  "% Vote for the Spanish Right, 1931",                  "spnxright_1931",
  "% Vote for the Catalan Left, 1931",                   "catxleft_1931",
  "% Vote for the Catalan Right, 1931",                  "catxright_1931",
  "Turnout in 1931",                                     "turnout_1931",
  
  # 1932 Catalan regional elections
  "% Vote for Catalan nationalist parties, 1932",                    "catxleft_1932 + catxright_1932",
  "% Vote for Left parties, 1932",                       "spnxleft_1932 + catxleft_1932",
  "% Vote for the Spanish Left, 1932",                   "spnxleft_1932",
  "% Vote for the Catalan Left, 1932",                   "catxleft_1932",
  "% Vote for the Catalan Right, 1932",                  "catxright_1932",
  "Turnout in 1932",                                     "turnout_1932",
  
  # 1933 Spanish national elections
  "% Vote for Catalan nationalist parties, 1933",                    "catxleft_1933 + catxright_1933",
  "% Vote for Left parties, 1933",                       "catxleft_1933 + spnxleft_1933",
  "% Vote for the Spanish Left, 1933",                   "spnxleft_1933",
  "% Vote for the Catalan Left, 1933",                   "catxleft_1933",
  "Turnout in 1933",                                     "turnout_1933",
  
  # February 1936 Spanish national elections
  "% Vote for Catalan nationalist parties, 1936",                    "catxleft_1936 + catxright_1936",
  "% Vote for Left parties, 1936",                       "catxleft_1936 + spnxleft_1936",
  "% Vote for the Spanish Left, 1936",                   "spnxleft_1936",
  "% Vote for the Catalan Left, 1936",                   "catxleft_1936",
  "% Vote for the Catalan Right, 1936",                  "catxright_1936",
  "Turnout in 1936",                                     "turnout_1936",
  
  # April 1936 Spanish national elections
  "% Vote for the Catalan Left, April 1936",             "catxleft_1936april",
  "% Vote for the Catalan Right, April 1936",            "catxright_1936april",
  "Turnout in April 1936",                               "turnout_1936april",
  
  # Membership in pre-war political associations
  "Catalanist political organizations, per 1000 people", "association_catalan_per1000",
  "Leftist political organizations, per 1000 people",    "association_left_per1000",
  
  # Membership in pre-war trade unions
  "Affiliated with the CNT, percent",                    "as.numeric(cntaffiliation)",
  "Affiliated with the UGT, percent",                    "as.numeric(ugtaffiliation)",
  
  # Population
  "Population in 1920",                                  "population_hecho_1920",
  "Population in 1930",                                  "population_hecho_1930",
  "Population change 1920 to 1930",                      "pop_dif_30_20",
  
  # Pre-front violence
  "Executions by the Left, per 1000 people",             "exec_by_left_prefront_per1000",
  "Church arson, per 1000 people",                       "burned_churches_per1000",
  "Church lootings 1936 to 1937, per 1000 people",       "looted_churches_36_37_per1000",
  "Other church damage, per 1000 people",                "undeter_damage_churches_per1000",
  "Priests murdered, per 1000 people",                   "priests_killed_per1000"
)

tab <- rows %>%
  rename(expr_col = expr, var_label = label) %>%
  mutate(out = map2(expr_col, var_label, ~ stat_row(.x, .y))) %>%  # pass safely
  unnest(out) %>%
  mutate(
    West = sprintf("%.1f (%.1f)", round(west_mean, 1), round(west_sd, 1)),
    East = sprintf("%.1f (%.1f)", round(east_mean, 1), round(east_sd, 1)),
    `t-statistic` = sprintf("%.2f", round(t_stat, 2))
  ) %>%
  select(var_label, West, East, `t-statistic`) %>%
  rename(label = var_label)

# Add N at the bottom for transparency
N <- aggregate_data %>% summarize(N = n()) %>% pull(N)


Table1 <- tab %>%
  kable(caption = "Prewar voting behavior, organizations, population, and pre-front violence, West vs. East",
        align = c("l","c","c","c"),
        col.names = c("", "West mean (sd)", "East mean (sd)", "t-statistic")) %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  add_header_above(c(" " = 1, " " = 3)) %>%
  footnote(general = paste0("N = ", N,
                            ". Two-sample Welch t-tests. Means and standard deviations shown by group."))

print(Table1)

# ------------------
# Joint balance regression and F-test
# ------------------
ortho_formula <- east_river ~ pop_dif_30_20 + population_hecho_1920 + population_hecho_1930 +
  exec_by_left_prefront_per1000 + burned_churches_per1000 + looted_churches_36_37_per1000 +
  undeter_damage_churches_per1000 + association_left_per1000 + association_catalan_per1000 +
  priests_killed_per1000 + as.numeric(cntaffiliation) + as.numeric(ugtaffiliation) +
  as.numeric(spnxleft_1931) + as.numeric(spnxright_1931) + as.numeric(catxleft_1931) + as.numeric(catxright_1931) +
  as.numeric(turnout_1931) +
  as.numeric(catxleft_1932) + as.numeric(catxright_1932) + as.numeric(spnxleft_1932) + as.numeric(turnout_1932) +
  as.numeric(spnxleft_1933) + as.numeric(catxleft_1933) + as.numeric(catxright_1933) +
  as.numeric(turnout_1933) +
  as.numeric(catxleft_1936) + as.numeric(catxright_1936) + as.numeric(spnxleft_1936) + as.numeric(turnout_1936) +
  as.numeric(catxleft_1936april) + as.numeric(catxright_1936april) + as.numeric(turnout_1936april)

ortho_fit <- lm(ortho_formula, data = aggregate_data)

# F-test of joint significance with p-value
fs   <- summary(ortho_fit)$fstatistic
Fval <- unname(fs["value"])
df1  <- unname(fs["numdf"])
df2  <- unname(fs["dendf"])
pval <- round(pf(Fval, df1, df2, lower.tail = FALSE),2)

cat(sprintf("F-test of joint significance: F(%d, %d) = %.2f, p-value = %.3f\n",
            df1, df2, Fval, pval))


# ------------------#
# ------------------#
#      Table 2      #
# ------------------#
# ------------------#

vars_upper <- c("massgrave_all", "massgrave_rep", "massgrave_nat")
labs_upper <- c("At least one mass grave",
                "At least one mass grave with Republicans' remains",
                "At least one mass grave with Nationalists' remains")

vars_lower <- c("bodycount_all", "bodycount_soldiers_rep", "bodycount_soldiers_nat")
labs_lower <- c("Number of bodies in mass graves",
                "Number of Republican bodies in mass graves",
                "Number of Nationalist bodies in mass graves")

digits <- c(
  massgrave_all = 2, massgrave_rep = 2, massgrave_nat = 2,   # proportions
  bodycount_all = 1, bodycount_soldiers_rep = 1, bodycount_soldiers_nat = 1  # counts
)

stat_compute <- function(var, label) {
  x <- aggregate_data %>% select(all_of(var), east_river)
  west <- x %>% filter(east_river == 0) %>% pull(1)
  east <- x %>% filter(east_river == 1) %>% pull(1)
  tt   <- t.test(west, east, alternative = "two.sided", var.equal = FALSE)
  stars <- dplyr::case_when(
    tt$p.value < 0.01 ~ "**",
    tt$p.value < 0.05 ~ "*",
    tt$p.value < 0.10 ~ "+",
    TRUE ~ ""
  )
  tibble::tibble(
    label = label,
    west_mean = mean(west, na.rm = TRUE),
    east_mean = mean(east, na.rm = TRUE),
    t_stat = unname(tt$statistic),
    stars = stars,
    var = var
  )
}

fmt_num <- function(x, d) sprintf(paste0("%.", d, "f"), round(x, d))

panel_upper <- purrr::map2_dfr(vars_upper, labs_upper, stat_compute) %>%
  dplyr::rename(var_name = var) %>%                    # avoid clash with base::var()
  dplyr::mutate(
    dig = digits[match(var_name, names(digits))],      # pick digits per variable
    West = fmt_num(west_mean, dig),
    East = fmt_num(east_mean, dig),
    `t-stat` = sprintf("%.2f%s", abs(round(t_stat, 2)), stars)
  ) %>%
  dplyr::select(label, West, East, `t-stat`)

panel_lower <- purrr::map2_dfr(vars_lower, labs_lower, stat_compute) %>%
  dplyr::rename(var_name = var) %>%                    # avoid clash with base::var()
  dplyr::mutate(
    dig = digits[match(var_name, names(digits))],      # pick digits per variable
    West = fmt_num(west_mean, dig),
    East = fmt_num(east_mean, dig),
    `t-stat` = sprintf("%.2f%s", abs(round(t_stat, 2)), stars)
  ) %>%
  dplyr::select(label, West, East, `t-stat`)

tbl <- bind_rows(
  tibble(label = "All violence:", West = "", East = "", `t-stat` = ""),
  panel_upper,
  tibble(label = "", West = "", East = "", `t-stat` = ""),
  tibble(label = "Number of bodies", West = "", East = "", `t-stat` = ""),
  panel_lower
)

N_val <- nrow(aggregate_data)

Table2 <- kable(tbl,
      format = "latex",
      booktabs = TRUE,
      align = c("l", "c", "c", "c"),
      col.names = c("", "West", "East", "t-stat"),
      caption = "Wartime Violence, Data on Mass Graves") %>%
  kable_styling(latex_options = c("hold_position")) %>%
  add_header_above(c(" " = 1, "All violence" = 3)) %>%
  footnote(
    general = paste0(
      "+ p<0.1; * p<0.05; ** p<0.01. P-values are based on two tailed t-tests ",
      "comparing east-west differences. N = ", N_val, "."
    ),
    threeparttable = TRUE,
    escape = FALSE
  )

print(Table2)

# ------------------#
# ------------------#
#   Text on p. 20   #
# ------------------#
# ------------------#

tt_bomb <- with(
  aggregate_data,
  t.test(
    bombardments_37_38[east_river == 0],
    bombardments_37_38[east_river == 1],
    alternative = "two.sided",
    var.equal = FALSE
  )
)

west_mean <- mean(aggregate_data$bombardments_37_38[aggregate_data$east_river == 0], na.rm = TRUE)
east_mean <- mean(aggregate_data$bombardments_37_38[aggregate_data$east_river == 1], na.rm = TRUE)
stars <- ifelse(tt_bomb$p.value < 0.01, "**",
                ifelse(tt_bomb$p.value < 0.05, "*",
                       ifelse(tt_bomb$p.value < 0.10, "+", "")))

sprintf(
  "West mean = %.2f, East mean = %.2f, t = %.2f%s, p = %.3f",
  west_mean, east_mean, unname(abs(tt_bomb$statistic)), stars, tt_bomb$p.value
)


# ------------------#
# ------------------#
#      Figure 2     #
# ------------------#
# ------------------#

# ------------------#
#      Panel A      #
# ------------------#

run_models <- function(years, prefix, suffix = NULL, data, variable_prefix = "model") {
  df <- data
  out <- vector("list", length(years))
  names(out) <- paste0(variable_prefix, "_", years)
  for (i in seq_along(years)) {
    y1 <- paste0(prefix, years[i])
    if (!is.null(suffix)) {
      y2 <- paste0(suffix, years[i])
      if (!(y1 %in% names(df)) || !(y2 %in% names(df))) {
        stop("Missing outcome columns for year ", years[i], ": ", y1, " or ", y2)
      }
      fml <- as.formula(paste0("I(", y1, " + ", y2, ") ~ east_river"))
    } else {
      if (!(y1 %in% names(df))) stop("Missing outcome column for year ", years[i], ": ", y1)
      fml <- as.formula(paste0(y1, " ~ east_river"))
    }
    dat_i <- df[, all.vars(fml), drop = FALSE]
    dat_i <- dat_i[complete.cases(dat_i), , drop = FALSE]
    out[[i]] <- lm(fml, data = dat_i)
  }
  out
}

extract_coefficients <- function(model_list, term_of_interest, model_label, year_labels) {
  coef_list <- vector("list", length(model_list))
  for (i in seq_along(model_list)) {
    coef_list[[i]] <- broom::tidy(model_list[[i]]) %>%
      dplyr::filter(term == term_of_interest) %>%
      dplyr::mutate(model = model_label, term = year_labels[i])
  }
  dplyr::bind_rows(coef_list)
}

congress_years <- c(1977, 1979, 1982, 1986, 1989, 1993, 1996, 2000, 2004, 2008, 2011, 2015, 2016, 2019)
cp_years       <- c(1980, 1984, 1988, 1992, 1995, 1999, 2003, 2006, 2010, 2012, 2015, 2017)

congress_years_labels <- paste0(congress_years, " - Congreso")
cp_years_labels       <- paste0(cp_years,       " - CP")

congress_models_left <- run_models(
  years = congress_years,
  prefix = "spnxleft_congress",
  suffix = "catxleft_congress",
  data = aggregate_data,
  variable_prefix = "left_congress"
)

cp_models_left <- run_models(
  years = cp_years,
  prefix = "spnxleft_cp",
  suffix = "catxleft_cp",
  data = aggregate_data,
  variable_prefix = "left_cp"
)

ols_left_congress_plot <- extract_coefficients(
  model_list = congress_models_left,
  term_of_interest = "east_river",
  model_label = "Left",
  year_labels = congress_years_labels
)

ols_left_cp_plot <- extract_coefficients(
  model_list = cp_models_left,
  term_of_interest = "east_river",
  model_label = "Left",
  year_labels = cp_years_labels
)

custom_order <- c(
  "1977 - Congreso",
  "1979 - Congreso",
  "1980 - CP",
  "1982 - Congreso",
  "1984 - CP",
  "1986 - Congreso",
  "1988 - CP",
  "1989 - Congreso",
  "1992 - CP",
  "1993 - Congreso",
  "1995 - CP",
  "1996 - Congreso",
  "1999 - CP",
  "2000 - Congreso",
  "2003 - CP",
  "2004 - Congreso",
  "2006 - CP",
  "2008 - Congreso",
  "2010 - CP",
  "2011 - Congreso",
  "2012 - CP",
  "2015 - Congreso",
  "2015 - CP",
  "2016 - Congreso",
  "2017 - CP",
  "2019 - Congreso"
)

ols_left_congress_plot <- extract_coefficients(
  model_list = congress_models_left,
  term_of_interest = "east_river",
  model_label = "Left",
  year_labels = congress_years_labels
)

ols_left_cp_plot <- extract_coefficients(
  model_list = cp_models_left,
  term_of_interest = "east_river",
  model_label = "Left",
  year_labels = cp_years_labels
)

all_left_plot <- bind_rows(ols_left_congress_plot, ols_left_cp_plot) %>%
  mutate(term = factor(term, levels = custom_order))

pd <- position_dodge(width = 0.2)

Figure2_PanelA <- ggplot(all_left_plot, aes(estimate, term)) +
  geom_point(size = 2, shape = 19, color = "steelblue", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, size = 0.5, color = "steelblue", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, size = 1.0, color = "steelblue", position = pd) +
  geom_vline(xintercept = 0, linetype = "dotted", color = "gray50", linewidth = 0.3) +
  theme_bw() +
  scale_x_continuous(
    "Estimated effect of living east of the river on voting for the Left",
    breaks = seq(-0.05, 0.20, by = 0.05),
    limits = c(-0.05, 0.16)
  ) +
  scale_y_discrete("") +
  coord_flip() +
  ggtitle("") +
  theme(
    legend.position = "none",
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
  )

print(Figure2_PanelA)


# ------------------#
#      Panel B      #
# ------------------#

catalan_congress_models <- run_models(
  years = congress_years,
  prefix = "catxleft_congress",
  suffix = "catxright_congress",
  data = aggregate_data,
  variable_prefix = "catalan_congress"
)

catalan_cp_models <- run_models(
  years = cp_years,
  prefix = "catxleft_cp",
  suffix = "catxright_cp",
  data = aggregate_data,
  variable_prefix = "catalan_cp"
)

ols_cat_cong <- extract_coefficients(
  model_list = catalan_congress_models,
  term_of_interest = "east_river",
  model_label = "Catalan (OLS)",
  year_labels = congress_years_labels
)

ols_cat_cp <- extract_coefficients(
  model_list = catalan_cp_models,
  term_of_interest = "east_river",
  model_label = "Catalan (OLS)",
  year_labels = cp_years_labels
)

cat_all <- bind_rows(ols_cat_cong, ols_cat_cp) %>%
  mutate(term = factor(term, levels = custom_order)) %>%
  arrange(term)

Figure2_PanelB <- ggplot(cat_all, aes(estimate, term)) +
  geom_point(size = 2, shape = 19, color = "#f58231", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, linewidth = 0.5, color = "#f58231", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, linewidth = 1.0, color = "#f58231", position = pd) +
  geom_vline(xintercept = 0, linetype = "dotted", color = "gray50", linewidth = 0.3) +
  scale_x_continuous("Estimated effect of living east of the river on voting for Catalan parties",
                     breaks = seq(-0.05, 0.20, by = 0.05),
                     limits = c(-0.05, 0.16),
                     labels = scales::number_format(accuracy = 0.01)) +
  scale_y_discrete("") +
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 90, vjust = 0.5),
        axis.text.y = element_text(angle = 0, size = 9)) +
  coord_flip()

print(Figure2_PanelB)

# ------------------#
#      Panel C      #
# ------------------#

catalan_left_congress_models <- run_models(
  years = congress_years,
  prefix = "catxleft_congress",
  data = aggregate_data,
  variable_prefix = "catalan_left_congress"
)

ols_catleft_congress_plot <- extract_coefficients(
  catalan_left_congress_models, 
  "east_river", 
  "Catalan Left", 
  congress_years_labels
)

catalan_left_cp_models <- run_models(
  years = cp_years,
  prefix = "catxleft_cp",
  data = aggregate_data,
  variable_prefix = "catalan_left_cp"
)

ols_catleft_cp_plot <- extract_coefficients(
  catalan_left_cp_models, "east_river", "Catalan Left", cp_years_labels
)

catalan_right_congress_models <- run_models(
  years = congress_years,
  prefix = "catxright_congress",
  data = aggregate_data,
  variable_prefix = "catalan_right_congress"
)

ols_catright_congress_plot <- extract_coefficients(
  catalan_right_congress_models, "east_river", "Catalan Right", congress_years_labels
)

catalan_right_cp_models <- run_models(
  years = cp_years,
  prefix = "catxright_cp",
  data = aggregate_data,
  variable_prefix = "catalan_right_cp"
)

ols_catright_cp_plot <- extract_coefficients(
  catalan_right_cp_models, "east_river", "Catalan Right", cp_years_labels
)

all_catleftright_plot <- bind_rows(
  ols_catleft_congress_plot,
  ols_catleft_cp_plot,
  ols_catright_congress_plot,
  ols_catright_cp_plot
) %>%
  mutate(
    term  = factor(term, levels = custom_order),
    model = factor(model, levels = c("Catalan Left","Catalan Right"))
  ) %>%
  arrange(term)

pd <- position_dodge(width = 0.25)

Figure2_PanelC <- ggplot(all_catleftright_plot, aes(estimate, term, color = model)) +
  geom_point(aes(shape = model), size = 2.2, position = pd, stroke = 1) +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, linewidth = 0.5, position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, linewidth = 1.0, position = pd) +
  geom_vline(xintercept = 0, linetype = "dotted", color = "gray50", linewidth = 0.3) +
  scale_color_manual(values = c("Catalan Left" = "darkgreen", "Catalan Right" = "darksalmon")) +
  scale_shape_manual(values = c("Catalan Left" = 17, "Catalan Right" = 19)) +  # ▲ vs ●
  scale_x_continuous(
    "Estimated effect of living east of the river on voting\nfor the Catalan Left and the Catalan Right",
    breaks = seq(-0.05, 0.20, by = 0.05),
    limits = c(-0.05, 0.16),
  ) +
  scale_y_discrete("") +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
    legend.background = element_rect(color = "grey70", fill = "white"),
    plot.title = element_text(face = "bold"),
    axis.text.x = element_text(angle = 90, vjust = 0.5),
    axis.text.y = element_text(angle = 0, size = 9)) +
  coord_flip()

print(Figure2_PanelC)


# ---------------------------------- #
# ---------------------------------- #
# ---------------------------------- #
#      SUPPLEMENTARY INFORMATION     #
# ---------------------------------- #
# ---------------------------------- #
# ---------------------------------- #

# --------------------------- #
# --------------------------- #
#  APPENDIX G: TABLE G.1      #
# --------------------------- #
# --------------------------- #

fmt_ms <- function(x) {
  m <- mean(x, na.rm = TRUE)
  s <- sd(x, na.rm = TRUE)
  sprintf("%.2f (%.2f)", m, s)
}

one_row <- function(df, var, nice_label) {
  west <- df %>% filter(east_river == 0) %>% pull({{ var }})
  east <- df %>% filter(east_river == 1) %>% pull({{ var }})
  
  tval <- tryCatch(
    t.test(west, east, alternative = "two.sided", var.equal = FALSE)$statistic[[1]],
    error = function(e) NA_real_
  )
  
  tibble(
    label = nice_label,
    west  = fmt_ms(west),
    east  = fmt_ms(east),
    t     = sprintf("%.2f", abs(tval))
  )
}

rows <- list(
  # Panel: All, frontline and postwar periods
  tibble(section = "All, frontline and postwar periods"),
  one_row(aggregate_data, prosec_after1938_exec_per1000, "Number of judicial proceedings initiated in 1938 or later resulting in Execution"),
  one_row(aggregate_data, prosec_after1938_impris_per1000, "Number of judicial proceedings initiated in 1938 or later resulting in Imprisonment"),
  one_row(aggregate_data, prosec_after1938_fine_per1000,     "Number of judicial proceedings initiated in 1938 or later resulting in Fine"),
  one_row(aggregate_data, prosec_after1938_profban_per1000,  "Number of judicial proceedings initiated in 1938 or later resulting in Professional ban"),

  # Panel: 1938, frontline period
  tibble(section = "1938, frontline period"),
  one_row(aggregate_data, prosec_1938_exec_per1000,     "Number of judicial proceedings initiated in 1938 resulting in Execution"),
  one_row(aggregate_data, prosec_1938_impris_per1000,  "Number of judicial proceedings initiated in 1938 resulting in Imprisonment"),
  one_row(aggregate_data, prosec_1938_fine_per1000,          "Number of judicial proceedings initiated in 1938 resulting in Fine"),
  one_row(aggregate_data, prosec_1938_profban_per1000,       "Number of judicial proceedings initiated in 1938 resulting in Professional ban"),

  # Panel: 1939, postwar period
  tibble(section = "1939, postwar period"),
  one_row(aggregate_data, prosec_1939_exec_per1000,     "Number of judicial proceedings initiated in 1939 resulting in Execution"),
  one_row(aggregate_data, prosec_1939_impris_per1000,  "Number of judicial proceedings initiated in 1939 resulting in Imprisonment"),
  one_row(aggregate_data, prosec_1939_fine_per1000,          "Number of judicial proceedings initiated in 1939 resulting in Fine"),
  one_row(aggregate_data, prosec_1939_profban_per1000,       "Number of judicial proceedings initiated in 1939 resulting in Professional ban"),

  # Panel: 1940 or later, postwar period
  tibble(section = "1940 or later, postwar period"),
  one_row(aggregate_data, prosec_1940_exec_per1000,     "Number of judicial proceedings initiated in 1940 or later resulting in Execution"),
  one_row(aggregate_data, prosec_1940_impris_per1000,  "Number of judicial proceedings initiated in 1940 or later resulting in Imprisonment"),
  one_row(aggregate_data, prosec_1940_fine_per1000,          "Number of judicial proceedings initiated in 1940 or later resulting in Fine"),
  one_row(aggregate_data, prosec_1940_profban_per1000,       "Number of judicial proceedings initiated in 1940 or later resulting in Professional ban")
)

table_g1 <- bind_rows(rows) %>%
  mutate(section = ifelse(is.na(west), section, NA_character_)) %>%
  rename(`West mean (sd)` = west,
         `East mean (sd)` = east,
         `t-statistic` = t)

print(table_g1, n = 24)


# --------------------------- #
# --------------------------- #
#  APPENDIX H.1: Figure H.1.  #
# --------------------------- #
# --------------------------- #

all_controls <- "exec_by_left_prefront_per1000 + burned_churches_per1000 + looted_churches_36_37_per1000 +
undeter_damage_churches_per1000 + priests_killed_per1000 + population_hecho_1920 + population_hecho_1930 +
association_left_per1000 + association_catalan_per1000 + as.numeric(cntaffiliation) + as.numeric(ugtaffiliation) +
as.numeric(spnxleft_1931) + as.numeric(spnxright_1931) + as.numeric(catxleft_1931) + as.numeric(catxright_1931) +
as.numeric(turnout_1931) + as.numeric(catxleft_1932) + as.numeric(catxright_1932) + as.numeric(spnxleft_1932) +
as.numeric(turnout_1932) + as.numeric(spnxleft_1933) + as.numeric(spnxright_1933) + as.numeric(catxleft_1933) +
as.numeric(catxright_1933) + as.numeric(turnout_1933) + as.numeric(catxleft_1936) + as.numeric(catxright_1936) +
as.numeric(spnxleft_1936) + as.numeric(turnout_1936)"


run_models_with_controls <- function(years, prefix, suffix = NULL, data, controls = NULL, variable_prefix = "model") {
    df <- data
    out <- vector("list", length(years))
    names(out) <- paste0(variable_prefix, "_", years)
    
    rhs <- if (is.null(controls) || nchar(trimws(controls)) == 0) {
      "east_river"
    } else {
      paste("east_river +", controls)
    }
    
    for (i in seq_along(years)) {
      y1 <- paste0(prefix, years[i])
      
      if (!is.null(suffix)) {
        y2 <- paste0(suffix, years[i])
        if (!(y1 %in% names(df)) || !(y2 %in% names(df))) stop("Missing outcome for ", years[i])
        
        outcome <- df[[y1]] + df[[y2]]
      } else {
        if (!(y1 %in% names(df))) stop("Missing outcome for ", years[i])
        outcome <- df[[y1]]
      }
      
      dat_i <- df[, c("east_river", all.vars(as.formula(paste("~", rhs)))), drop = FALSE]
      dat_i$outcome <- outcome
      
      dat_i <- stats::na.omit(dat_i)
      
      fml <- stats::as.formula(paste("outcome ~", rhs))
      out[[i]] <- stats::lm(fml, data = dat_i)
    }
    out
  }

# ------------------#
#      Panel A      #
# ------------------#

congress_models_left_controls <- run_models_with_controls(
  years = congress_years,
  prefix = "spnxleft_congress",
  suffix = "catxleft_congress",
  data   = aggregate_data,
  controls = all_controls,
  variable_prefix = "left_congress"
)

cp_models_left_controls <- run_models_with_controls(
  years = cp_years,
  prefix = "spnxleft_cp",
  suffix = "catxleft_cp",
  data   = aggregate_data,
  controls = all_controls,
  variable_prefix = "left_cp"
)

ols_left_congress_controls <- extract_coefficients(
  model_list     = congress_models_left_controls,
  term_of_interest = "east_river",
  model_label    = "Left",
  year_labels    = congress_years_labels
)

ols_left_cp_controls <- extract_coefficients(
  model_list     = cp_models_left_controls,
  term_of_interest = "east_river",
  model_label    = "Left",
  year_labels    = cp_years_labels
)

all_left_plot_controls <- bind_rows(ols_left_congress_controls, ols_left_cp_controls) %>%
  mutate(term = factor(term, levels = custom_order))

pd <- position_dodge(width = 0.2)

FigureH1_PanelA <- ggplot(all_left_plot_controls, aes(estimate, term)) +
  geom_point(size = 2, shape = 19, color = "steelblue", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, size = 0.5, color = "steelblue", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, size = 1.0, color = "steelblue", position = pd) +
  geom_vline(xintercept = 0, linetype = "dotted", color = "gray50", linewidth = 0.3) +
  theme_bw() +
  scale_x_continuous(
    "Estimated effect of living east of the river on voting for the Left (with controls)",
    breaks = seq(-0.05, 0.20, by = 0.05),
    limits = c(-0.05, 0.16)
  ) +
  scale_y_discrete("") +
  coord_flip() +
  ggtitle("") +
  theme(
    legend.position = "none",
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
  )

print(FigureH1_PanelA)

# ------------------#
#      Panel B      #
# ------------------#

catalan_congress_models_controls <- run_models_with_controls(
  years = congress_years,
  prefix = "catxleft_congress",
  suffix = "catxright_congress",
  data   = aggregate_data,
  controls = all_controls,
  variable_prefix = "catalan_congress"
)

catalan_cp_models_controls <- run_models_with_controls(
  years = cp_years,
  prefix = "catxleft_cp",
  suffix = "catxright_cp",
  data   = aggregate_data,
  controls = all_controls,
  variable_prefix = "catalan_cp"
)

ols_cat_cong_controls <- extract_coefficients(
  model_list     = catalan_congress_models_controls,
  term_of_interest = "east_river",
  model_label    = "Catalan (OLS)",
  year_labels    = congress_years_labels
)

ols_cat_cp_controls <- extract_coefficients(
  model_list     = catalan_cp_models_controls,
  term_of_interest = "east_river",
  model_label    = "Catalan (OLS)",
  year_labels    = cp_years_labels
)

cat_all_controls <- bind_rows(ols_cat_cong_controls, ols_cat_cp_controls) %>%
  mutate(term = factor(term, levels = custom_order)) %>%
  arrange(term)

FigureH1_PanelB <- ggplot(cat_all_controls, aes(estimate, term)) +
  geom_point(size = 2, shape = 19, color = "#f58231", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, linewidth = 0.5, color = "#f58231", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, linewidth = 1.0, color = "#f58231", position = pd) +
  geom_vline(xintercept = 0, linetype = "dotted", color = "gray50", linewidth = 0.3) +
  scale_x_continuous("Estimated effect of living east of the river on voting for Catalan parties",
                     breaks = seq(-0.05, 0.20, by = 0.05),
                     limits = c(-0.05, 0.16),
                     labels = scales::number_format(accuracy = 0.01)) +
  scale_y_discrete("") +
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 90, vjust = 0.5),
        axis.text.y = element_text(angle = 0, size = 9)) +
  coord_flip()

print(FigureH1_PanelB)

# --------------------------- #
# --------------------------- #
#          APPENDIX I         #
# --------------------------- #
# --------------------------- #

make_block <- function(data, years, lab, left_sp_prefix, left_cat_prefix, right_sp_prefix, right_cat_prefix) {
  purrr::map_dfr(years, function(y) {
    tibble(
      Year = sprintf("%d - %s", y, lab),
      `Spanish Left`  = mean(data[[paste0(left_sp_prefix,  y)]],  na.rm = TRUE),
      `Catalan Left`  = mean(data[[paste0(left_cat_prefix, y)]],  na.rm = TRUE),
      `Spanish Right` = mean(data[[paste0(right_sp_prefix, y)]],  na.rm = TRUE),
      `Catalan Right` = mean(data[[paste0(right_cat_prefix, y)]], na.rm = TRUE)
    )
  })
}

tab_congress <- make_block(
  data = aggregate_data, years = congress_years, lab = "Congreso",
  left_sp_prefix = "spnxleft_congress",
  left_cat_prefix = "catxleft_congress",
  right_sp_prefix = "spnxright_congress",
  right_cat_prefix = "catxright_congress"
)

tab_cp <- make_block(
  data = aggregate_data, years = cp_years, lab = "Catalan Parliament",
  left_sp_prefix = "spnxleft_cp",
  left_cat_prefix = "catxleft_cp",
  right_sp_prefix = "spnxright_cp",
  right_cat_prefix = "catxright_cp"
)

fmt_pct <- function(x) sprintf("%.1f%%", 100 * x)

tab_congress_fmt <- tab_congress %>%
  mutate(across(-Year, fmt_pct))

tab_cp_fmt <- tab_cp %>%
  mutate(across(-Year, fmt_pct))

print(tab_congress_fmt)
print(tab_cp_fmt)

# --------------------------- #
# --------------------------- #
#   APPENDIX J: Figure J.1    #
# --------------------------- #
# --------------------------- #

ols_catleft_congress_plot <- extract_coefficients(
  model_list = catalan_left_congress_models,
  term_of_interest = "east_river",
  model_label = "Catalan Left",
  year_labels = congress_years_labels
)

spanish_left_congress_models <- run_models(
  years = congress_years,
  prefix = "spnxleft_congress",
  data = aggregate_data,
  variable_prefix = "spanish_left_congress"
)

ols_spaleft_congress_plot <- extract_coefficients(
  model_list = spanish_left_congress_models,
  term_of_interest = "east_river",
  model_label = "Spanish Left",
  year_labels = congress_years_labels
)

ols_catleft_cp_plot <- extract_coefficients(
  model_list = catalan_left_cp_models,
  term_of_interest = "east_river",
  model_label = "Catalan Left",
  year_labels = cp_years_labels
)

spanish_left_cp_models <- run_models(
  years = cp_years,
  prefix = "spnxleft_cp",
  data = aggregate_data,
  variable_prefix = "spanish_left_cp"
)

ols_spaleft_cp_plot <- extract_coefficients(
  model_list = spanish_left_cp_models,
  term_of_interest = "east_river",
  model_label = "Spanish Left",
  year_labels = cp_years_labels
)

custom_order <- c(
  "1977 - Congreso","1979 - Congreso","1980 - CP","1982 - Congreso","1984 - CP",
  "1986 - Congreso","1988 - CP","1989 - Congreso","1992 - CP","1993 - Congreso",
  "1995 - CP","1996 - Congreso","1999 - CP","2000 - Congreso","2003 - CP",
  "2004 - Congreso","2006 - CP","2008 - Congreso","2010 - CP","2011 - Congreso",
  "2012 - CP","2015 - Congreso","2015 - CP","2016 - Congreso","2017 - CP",
  "2019 - Congreso"
)

all_left_catspa_plot <- bind_rows(
  ols_catleft_congress_plot,
  ols_spaleft_congress_plot,
  ols_catleft_cp_plot,
  ols_spaleft_cp_plot
) %>%
  mutate(
    term  = factor(term, levels = custom_order),
    model = factor(model, levels = c("Catalan Left","Spanish Left"))
  ) %>%
  arrange(term)

pd <- position_dodge(width = 0.35)

FigureJ <- ggplot(all_left_catspa_plot,
                                               aes(x = term, y = estimate,
                                                   color = model, shape = model)) +
  geom_point(size = 2.2, position = pd) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error),
                width = 0, linewidth = 0.5, position = pd) +
  geom_errorbar(aes(ymin = estimate - 1.645 * std.error,
                    ymax = estimate + 1.645 * std.error),
                width = 0, linewidth = 1.0, position = pd) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "gray50", linewidth = 0.3) +
  scale_color_manual(values = c("Catalan Left" = "darkgreen",
                                "Spanish Left" = "red")) +
  scale_shape_manual(values = c("Catalan Left" = 17, "Spanish Left" = 19)) +
  scale_y_continuous("Estimated effect of living east of the river on voting for \nCatalan Left and Spanish Left parties",
                     breaks = seq(-0.10, 0.10, by = 0.05),
                     limits = c(-0.10, 0.10)) +
  scale_x_discrete("") +
  theme_bw() +
  theme(
    legend.position = "top",
    legend.title = element_blank(),
    legend.background = element_rect(colour = "grey50", fill = "white"),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
  )

print(FigureJ)

# --------------------------- #
# --------------------------- #
#         APPENDIX K          #
# --------------------------- #
# --------------------------- #

starify <- function(p) {
  case_when(
    is.na(p)        ~ "",
    p < .01         ~ "**",
    p < .05         ~ "*",
    p < .10         ~ "\u207A",  # superscript plus
    TRUE            ~ ""
  )
}

fmt_cell <- function(est, se, p) {
  paste0(sprintf("%.3f", est), starify(p), " (", sprintf("%.3f", se), ")")
}

make_table_df <- function(model_list, years_vec) {
  
  td <- map2_dfr(model_list, years_vec, function(m, yr) {
    t <- tidy(m) %>% select(term, estimate, std.error, p.value)
    n <- glance(m)$nobs[1]
    t %>% mutate(year = yr, nobs = n)
  }) %>% arrange(year)
  
  yr_df <- tibble(year = years_vec)
  
  east <- td %>% filter(term == "east_river") %>%
    transmute(year, cell = fmt_cell(estimate, std.error, p.value)) %>%
    right_join(yr_df, by = "year") %>% mutate(cell = ifelse(is.na(cell), "", cell))
  
  cons <- td %>% filter(term == "(Intercept)") %>%
    transmute(year, cell = fmt_cell(estimate, std.error, p.value)) %>%
    right_join(yr_df, by = "year") %>% mutate(cell = ifelse(is.na(cell), "", cell))
  
  obs  <- td %>% distinct(year, nobs) %>%
    right_join(yr_df, by = "year") %>% mutate(nobs = ifelse(is.na(nobs), NA_integer_, nobs))
  
  yrs <- as.character(years_vec)
  
  east_row  <- data.frame(Row = "East of the river", t(east$cell), check.names = FALSE)
  names(east_row)[-1] <- yrs
  
  cons_row  <- data.frame(Row = "Constant", t(cons$cell), check.names = FALSE)
  names(cons_row)[-1] <- yrs
  
  obs_row   <- data.frame(Row = "Observations", t(as.character(obs$nobs)), check.names = FALSE)
  names(obs_row)[-1] <- yrs
  
  out <- bind_rows(east_row, cons_row, obs_row)
  
  as.data.frame(out)
}

# --------------------------- #
#         TABLE K.1           #
# --------------------------- #

TableK1 <- make_table_df(congress_models_left, congress_years)

print(TableK1, row.names = FALSE)

# --------------------------- #
#         TABLE K.2           #
# --------------------------- #

TableK2 <- make_table_df(cp_models_left, cp_years)

print(TableK2, row.names = FALSE)

# --------------------------- #
#         TABLE K.3           #
# --------------------------- #

TableK3 <- make_table_df(catalan_congress_models, congress_years)

print(TableK3, row.names = FALSE)

# --------------------------- #
#         TABLE K.4           #
# --------------------------- #

TableK4 <- make_table_df(catalan_cp_models, cp_years)

print(TableK4, row.names = FALSE)

# --------------------------- #
#         TABLE K.5           #
# --------------------------- #

TableK5 <- make_table_df(catalan_left_congress_models, congress_years)

print(TableK5, row.names = FALSE)

# --------------------------- #
#         TABLE K.6           #
# --------------------------- #

TableK6 <- make_table_df(catalan_left_cp_models, cp_years)

print(TableK6, row.names = FALSE)

# --------------------------- #
#         TABLE K.7           #
# --------------------------- #

TableK7 <- make_table_df(spanish_left_congress_models, congress_years)

print(TableK7, row.names = FALSE)

# --------------------------- #
#         TABLE K.8           #
# --------------------------- #

TableK8 <- make_table_df(spanish_left_cp_models, cp_years)

print(TableK8, row.names = FALSE)

# --------------------------- #
#         TABLE K.9           #
# --------------------------- #

TableK9 <- make_table_df(catalan_right_congress_models, congress_years)

print(TableK9, row.names = FALSE)

# --------------------------- #
#        TABLE K.10           #
# --------------------------- #

TableK10 <- make_table_df(catalan_right_cp_models, cp_years)

print(TableK10, row.names = FALSE)

# --------------------------- #
# --------------------------- #
#         APPENDIX L          #
# --------------------------- #
# --------------------------- #

run_iv_models <-
  function(years, prefix, suffix = NULL, variable_prefix, instrument, instrumented, control) {
    models <- list()
    
    first_stage_formula <- as.formula(
      paste0(instrumented, " ~ ", instrument, " + ", control)
    )
    
    for (year in years) {
      if (!is.null(suffix)) {
        outcome_formula <- as.formula(
          paste0(
            "I(", prefix, year, " + ", suffix, year, ") ~ ",
            instrumented, " + ", control,
            " | . - ", instrumented, " + ", instrument
          )
        )
      } else {
        outcome_formula <- as.formula(
          paste0(
            prefix, year, " ~ ",
            instrumented, " + ", control,
            " | . - ", instrumented, " + ", instrument
          )
        )
      }
      
      iv_fit <- ivreg(outcome_formula, data = aggregate_data)
      
      fs_fit <- lm(first_stage_formula, data = aggregate_data)
      
      fs_drop1 <- tryCatch(drop1(fs_fit, test = "F"), error = function(e) NULL)
      fs_F_value <- NA_real_
      fs_p_value <- NA_real_
      if (!is.null(fs_drop1) && instrument %in% rownames(fs_drop1)) {
        fs_F_value <- fs_drop1[instrument, "F value"]
        fs_p_value <- fs_drop1[instrument, "Pr(>F)"]
      }
      
      model_name <- paste0("ivmod_", variable_prefix, year)
      models[[model_name]] <- list(
        iv = iv_fit,
        first_stage = fs_fit,
        first_stage_F = fs_F_value,
        first_stage_p = fs_p_value
      )
      
      summary(iv_fit)
      summary(fs_fit)
      if (!is.na(fs_F_value)) {
        message(
          sprintf(
            "[%s] First-stage F on %s: %.3f, p-value: %.3g",
            model_name, instrument, fs_F_value, fs_p_value
          )
        )
      }
    }
    
    return(models)
  }

unwrap_models <- function(model_list) {
  lapply(model_list, function(m) if (inherits(m, "ivreg")) m else m$iv)
}

# --------------------------- #
#        FIGURE L.1           #
# --------------------------- #

congress_iv_models_left <- run_iv_models(
  years = congress_years,
  prefix = "spnxleft_congress",
  suffix = "catxleft_congress",
  variable_prefix = "left_congress",
  instrument = "east_river",
  instrumented = "violence_nat",
  control = "violence_rep"
)

cp_iv_models_left <- run_iv_models(
  years = cp_years,
  prefix = "spnxleft_cp",
  suffix = "catxleft_cp",
  variable_prefix = "left_cp",
  instrument = "east_river",
  instrumented = "violence_nat",
  control = "violence_rep"
)

catalan_congress_iv_models <- run_iv_models(
  years = congress_years,
  prefix = "catxleft_congress",
  suffix = "catxright_congress",
  variable_prefix = "catalan_congress",
  instrument = "east_river",
  instrumented = "violence_nat",
  control = "violence_rep"
)

catalan_cp_iv_models <- run_iv_models(
  years = cp_years,
  prefix = "catxleft_cp",
  suffix = "catxright_cp",
  variable_prefix = "catalan_cp",
  instrument = "east_river",
  instrumented = "violence_nat",
  control = "violence_rep"
)

iv_left_congress_plot <- extract_coefficients(
  model_list = unwrap_models(congress_iv_models_left),
  term_of_interest = "violence_nat",
  model_label = "Left (IV)",
  year_labels = congress_years_labels
)

iv_left_cp_plot <- extract_coefficients(
  model_list = unwrap_models(cp_iv_models_left),
  term_of_interest = "violence_nat",
  model_label = "Left (IV)",
  year_labels = cp_years_labels
)

iv_cat_congress_plot <- extract_coefficients(
  model_list = unwrap_models(catalan_congress_iv_models),
  term_of_interest = "violence_nat",
  model_label = "Catalan",
  year_labels = congress_years_labels
)

iv_cat_cp_plot <- extract_coefficients(
  model_list = unwrap_models(catalan_cp_iv_models),
  term_of_interest = "violence_nat",
  model_label = "Catalan",
  year_labels = cp_years_labels
)

all_left_iv_plot <- bind_rows(iv_left_congress_plot, iv_left_cp_plot) %>%
  arrange(desc(term))

all_left_iv_plot <- all_left_iv_plot %>%
  arrange(desc(term))

all_cat_iv_plot <- bind_rows(
  iv_cat_congress_plot,
  iv_cat_cp_plot
)

all_cat_iv_plot <- all_cat_iv_plot %>%
  arrange(desc(term))

pd <- position_dodge(width = 0.2)

FigureL1a <- ggplot(all_left_iv_plot, aes(estimate, term, color = model)) +
  geom_point(aes(shape = model), size = 2, position = pd) +
  scale_color_manual(name = "Model", values = c("steelblue")) +
  scale_shape_manual(name = "Model", values = c(17)) +
  theme_bw() +
  scale_x_continuous("Estimated effect of Nationalist violence on voting for the Left",
                     breaks = seq(-0.5, 0.5, by = 0.1), labels = seq(-0.5, 0.5, by = 0.1)) +
  scale_y_discrete("") +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, size = 0.5, position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, size = 1, position = pd) +
  geom_vline(xintercept = 0, linetype = "dotted",
             color = "gray50", linewidth = 0.3) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  coord_flip()

FigureL1b <- ggplot(all_cat_iv_plot, aes(estimate, term, color = model)) +
  geom_point(aes(shape = model), size = 3, position = pd) +
  scale_color_manual(name = "Model", values = c("coral")) +
  scale_shape_manual(name = "Model", values = c(19)) +
  theme_bw() +
  scale_x_continuous("Estimated effect of Nationalist violence on voting for Catalan parties",
                     breaks = seq(-0.5, 0.8, by = 0.1), labels = seq(-0.5, 0.8, by = 0.1)) +
  scale_y_discrete("") +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, size = 0.5, position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, size = 1, position = pd) +
  geom_vline(xintercept = 0, linetype = "dotted",
             color = "gray50", size = 0.5) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  coord_flip()


print(FigureL1a)
print(FigureL1b)

# --------------------------- #
#        FIGURE L.2           #
# --------------------------- #

congress_iv_models_left_civ <- run_iv_models(
  years = congress_years,
  prefix = "spnxleft_congress",
  suffix = "catxleft_congress",
  variable_prefix = "left_congress_civ",
  instrument = "east_river",
  instrumented = "violence_nat_civ",
  control = "violence_rep_civ"
)

cp_iv_models_left_civ <- run_iv_models(
  years = cp_years,
  prefix = "spnxleft_cp",
  suffix = "catxleft_cp",
  variable_prefix = "left_cp_civ",
  instrument = "east_river",
  instrumented = "violence_nat_civ",
  control = "violence_rep_civ"
)

catalan_congress_iv_models_civ <- run_iv_models(
  years = congress_years,
  prefix = "catxleft_congress",
  suffix = "catxright_congress",
  variable_prefix = "catalan_congress_civ",
  instrument = "east_river",
  instrumented = "violence_nat_civ",
  control = "violence_rep_civ"
)

catalan_cp_iv_models_civ <- run_iv_models(
  years = cp_years,
  prefix = "catxleft_cp",
  suffix = "catxright_cp",
  variable_prefix = "catalan_cp_civ",
  instrument = "east_river",
  instrumented = "violence_nat_civ",
  control = "violence_rep_civ"
)

iv_left_congress_plot_civ <- extract_coefficients(
  model_list = unwrap_models(congress_iv_models_left_civ),
  term_of_interest = "violence_nat_civ",
  model_label = "Left, IV",
  year_labels = congress_years_labels
)

iv_left_cp_plot_civ <- extract_coefficients(
  model_list = unwrap_models(cp_iv_models_left_civ),
  term_of_interest = "violence_nat_civ",
  model_label = "Left, IV",
  year_labels = cp_years_labels
)

iv_cat_congress_plot_civ <- extract_coefficients(
  model_list = unwrap_models(catalan_congress_iv_models_civ),
  term_of_interest = "violence_nat_civ",
  model_label = "Catalan",
  year_labels = congress_years_labels
)

iv_cat_cp_plot_civ <- extract_coefficients(
  model_list = unwrap_models(catalan_cp_iv_models_civ),
  term_of_interest = "violence_nat_civ",
  model_label = "Catalan",
  year_labels = cp_years_labels
)

all_left_iv_plot_civ <- bind_rows(iv_left_congress_plot_civ, iv_left_cp_plot_civ) %>%
  arrange(desc(term))

all_cat_iv_plot_civ <- bind_rows(iv_cat_congress_plot_civ, iv_cat_cp_plot_civ) %>%
  arrange(desc(term))

pd <- position_dodge(width = 0.2)

FigureL3a <- ggplot(all_left_iv_plot_civ, aes(estimate, term, color = model)) +
  geom_point(aes(shape = model), size = 2, position = pd) +
  scale_color_manual(name = "Model", values = c("steelblue")) +
  scale_shape_manual(name = "Model", values = c(17)) +
  theme_bw() +
  scale_x_continuous("Estimated effect of Nationalist civilian violence on voting for the Left",
                     breaks = seq(-0.5, 0.5, by = 0.1), labels = seq(-0.5, 0.5, by = 0.1)) +
  scale_y_discrete("") +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, size = 0.5, position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, size = 1, position = pd) +
  geom_vline(xintercept = 0, linetype = "dotted",
             color = "gray50", linewidth = 0.3) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  coord_flip()

FigureL3b <- ggplot(all_cat_iv_plot_civ, aes(estimate, term, color = model)) +
  geom_point(aes(shape = model), size = 3, position = pd) +
  scale_color_manual(name = "Model", values = c("coral")) +
  scale_shape_manual(name = "Model", values = c(19)) +
  theme_bw() +
  scale_x_continuous("Estimated effect of Nationalist civilian violence on voting for Catalan parties",
                     breaks = seq(-0.5, 0.8, by = 0.1), labels = seq(-0.5, 0.8, by = 0.1)) +
  scale_y_discrete("") +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, size = 0.5, position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, size = 1, position = pd) +
  geom_vline(xintercept = 0, linetype = "dotted",
             color = "gray50", size = 0.5) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  coord_flip()

print(FigureL3a)
print(FigureL3b)

# --------------------------- #
# --------------------------- #
#         APPENDIX O          #
# --------------------------- #
# --------------------------- #

lower_lat <- 40.9
upper_lat <- 41.3

aggregate_data_exclEbro <- aggregate_data[
  which(as.numeric(aggregate_data$centroid_lat) < lower_lat |
          as.numeric(aggregate_data$centroid_lat) > upper_lat),
]

# --------------------------- #
#         TABLE O.1           #
# --------------------------- #

vars_upper <- c("massgrave_all", "massgrave_rep", "massgrave_nat")
labs_upper <- c(
  "At least one mass grave",
  "At least one mass grave with Republicans' remains",
  "At least one mass grave with Nationalists' remains"
)

vars_lower <- c("bodycount_all", "bodycount_soldiers_rep", "bodycount_soldiers_nat")
labs_lower <- c(
  "Number of bodies in mass graves",
  "Number of Republican bodies in mass graves",
  "Number of Nationalist bodies in mass graves"
)

digits <- c(
  massgrave_all = 2, massgrave_rep = 2, massgrave_nat = 2,
  bodycount_all = 1, bodycount_soldiers_rep = 1, bodycount_soldiers_nat = 1
)

stat_compute_excl <- function(var, label) {
  x <- aggregate_data_exclEbro %>% dplyr::select(dplyr::all_of(var), east_river)
  west <- x %>% dplyr::filter(east_river == 0) %>% dplyr::pull(1)
  east <- x %>% dplyr::filter(east_river == 1) %>% dplyr::pull(1)
  tt   <- t.test(west, east, alternative = "two.sided", var.equal = FALSE)
  stars <- dplyr::case_when(
    tt$p.value < 0.01 ~ "**",
    tt$p.value < 0.05 ~ "*",
    tt$p.value < 0.10 ~ "+",
    TRUE ~ ""
  )
  tibble::tibble(
    label = label,
    west_mean = mean(west, na.rm = TRUE),
    east_mean = mean(east, na.rm = TRUE),
    t_stat = unname(tt$statistic),
    stars = stars,
    var = var
  )
}

fmt_num <- function(x, d) sprintf(paste0("%.", d, "f"), round(x, d))

panel_upper_excl <- purrr::map2_dfr(vars_upper, labs_upper, stat_compute_excl) %>%
  dplyr::rename(var_name = var) %>%
  dplyr::mutate(
    dig = digits[match(var_name, names(digits))],
    West = fmt_num(west_mean, dig),
    East = fmt_num(east_mean, dig),
    `t-stat` = sprintf("%.2f%s", abs(round(t_stat, 2)), stars)
  ) %>%
  dplyr::select(label, West, East, `t-stat`)

panel_lower_excl <- purrr::map2_dfr(vars_lower, labs_lower, stat_compute_excl) %>%
  dplyr::rename(var_name = var) %>%
  dplyr::mutate(
    dig = digits[match(var_name, names(digits))],
    West = fmt_num(west_mean, dig),
    East = fmt_num(east_mean, dig),
    `t-stat` = sprintf("%.2f%s", abs(round(t_stat, 2)), stars)
  ) %>%
  dplyr::select(label, West, East, `t-stat`)

tbl_O1_excl <- dplyr::bind_rows(
  tibble::tibble(label = "All violence:", West = "", East = "", `t-stat` = ""),
  panel_upper_excl,
  tibble::tibble(label = "", West = "", East = "", `t-stat` = ""),
  tibble::tibble(label = "Number of bodies in mass graves", West = "", East = "", `t-stat` = ""),
  panel_lower_excl
)

N_val_excl <- nrow(aggregate_data_exclEbro)

TableO1 <- knitr::kable(
  tbl_O1_excl,
  format = "latex",
  booktabs = TRUE,
  align = c("l", "c", "c", "c"),
  col.names = c("", "West", "East", "t-stat"),
  caption = "Appendix O.1, objective measures of wartime experiences across the river, excluding towns affected by the Republican Ebro offensive"
) %>%
  kableExtra::kable_styling(latex_options = c("hold_position")) %>%
  kableExtra::add_header_above(c(" " = 1, "All violence" = 3)) %>%
  kableExtra::footnote(
    general = paste0(
      "+ p<0.1; * p<0.05; ** p<0.01. P-values are based on two tailed t-tests ",
      "comparing east, west differences. N = ", N_val_excl, ". ",
      "The sample excludes municipalities with centroid latitude between ",
      lower_lat, " and ", upper_lat, "."
    ),
    threeparttable = TRUE,
    escape = FALSE
  )

print(TableO1)

# --------------------------- #
#        FIGURE O.3           #
# --------------------------- #

# --------------------------- #
#          PANEL A            #
# --------------------------- #

congress_models_left_exclEbro <- run_models(
  years = congress_years,
  prefix = "spnxleft_congress",
  suffix = "catxleft_congress",
  data = aggregate_data_exclEbro,
  variable_prefix = "left_congress_exclEbro"
)

cp_models_left_exclEbro <- run_models(
  years = cp_years,
  prefix = "spnxleft_cp",
  suffix = "catxleft_cp",
  data = aggregate_data_exclEbro,
  variable_prefix = "left_cp_exclEbro"
)

ols_left_congress_plot_exclEbro <- extract_coefficients(
  model_list = congress_models_left_exclEbro,
  term_of_interest = "east_river",
  model_label = "Left, excl. Ebro",
  year_labels = congress_years_labels
)

ols_left_cp_plot_exclEbro <- extract_coefficients(
  model_list = cp_models_left_exclEbro,
  term_of_interest = "east_river",
  model_label = "Left, excl. Ebro",
  year_labels = cp_years_labels
)

all_left_plot_exclEbro <- bind_rows(
  ols_left_congress_plot_exclEbro,
  ols_left_cp_plot_exclEbro
) %>%
  mutate(term = factor(term, levels = custom_order))

pd_excl <- position_dodge(width = 0.2)

FigureO1_PanelA <- ggplot(all_left_plot_exclEbro, aes(estimate, term)) +
  geom_point(size = 2, shape = 19, color = "steelblue", position = pd_excl) +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, size = 0.5, color = "steelblue", position = pd_excl) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, size = 1.0, color = "steelblue", position = pd_excl) +
  geom_vline(xintercept = 0, linetype = "dotted", color = "gray50", linewidth = 0.3) +
  theme_bw() +
  scale_x_continuous(
    "Estimated effect of living east of the river on voting for the Left,\n excluding towns affected by the Ebro offensive",
    breaks = seq(-0.05, 0.20, by = 0.05),
    limits = c(-0.07, 0.16)
  ) +
  scale_y_discrete("") +
  coord_flip() +
  ggtitle("") +
  theme(
    legend.position = "none",
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
  )

print(FigureO1_PanelA)

# --------------------------- #
#          PANEL B            #
# --------------------------- #

catalan_congress_models_exclEbro <- run_models(
  years = congress_years,
  prefix = "catxleft_congress",
  suffix = "catxright_congress",
  data = aggregate_data_exclEbro,
  variable_prefix = "catalan_congress_exclEbro"
)

catalan_cp_models_exclEbro <- run_models(
  years = cp_years,
  prefix = "catxleft_cp",
  suffix = "catxright_cp",
  data = aggregate_data_exclEbro,
  variable_prefix = "catalan_cp_exclEbro"
)

ols_cat_cong_exclEbro <- extract_coefficients(
  model_list = catalan_congress_models_exclEbro,
  term_of_interest = "east_river",
  model_label = "Catalan, excl. Ebro",
  year_labels = congress_years_labels
)

ols_cat_cp_exclEbro <- extract_coefficients(
  model_list = catalan_cp_models_exclEbro,
  term_of_interest = "east_river",
  model_label = "Catalan, excl. Ebro",
  year_labels = cp_years_labels
)

cat_all_exclEbro <- bind_rows(ols_cat_cong_exclEbro, ols_cat_cp_exclEbro) %>%
  mutate(term = factor(term, levels = custom_order)) %>%
  arrange(term)

FigureO1_PanelB <- ggplot(cat_all_exclEbro, aes(estimate, term)) +
  geom_point(size = 2, shape = 19, color = "#f58231", position = pd_excl) +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, linewidth = 0.5, color = "#f58231", position = pd_excl) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, linewidth = 1.0, color = "#f58231", position = pd_excl) +
  geom_vline(xintercept = 0, linetype = "dotted", color = "gray50", linewidth = 0.3) +
  scale_x_continuous(
    "Estimated effect of living east of the river on voting for Catalan parties,\n excluding towns affected by the Ebro offensive",
    breaks = seq(-0.05, 0.20, by = 0.05),
    limits = c(-0.06, 0.16),
    labels = scales::number_format(accuracy = 0.01)
  ) +
  scale_y_discrete("") +
  theme_bw() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 90, vjust = 0.5),
    axis.text.y = element_text(angle = 0, size = 9)
  ) +
  coord_flip()

print(FigureO1_PanelB)

# --------------------------- #
# --------------------------- #
#         APPENDIX P          #
# --------------------------- #
# --------------------------- #

# --------------------------------- #
# TABLE P.1: Panel "Municipalities" #
# --------------------------------- #

bandwidths <- c(20, 17, 12, 7)

muni_counts <- purrr::map_dfr(bandwidths, function(bw) {
  d <- aggregate_data %>%
    dplyr::filter(!is.na(dist_river_combined),
                  dist_river_combined <= bw,
                  !is.na(east_river)) %>%
    dplyr::mutate(east_flag = ifelse(east_river %in% c(1, "1", TRUE, "TRUE"), 1L, 0L))
  
  west_n <- sum(d$east_flag == 0L, na.rm = TRUE)
  east_n <- sum(d$east_flag == 1L, na.rm = TRUE)
  
  tibble::tibble(
    `Bandwidth (km)`     = bw,
    `West of the Line`   = west_n,
    `East of the Line`   = east_n,
    `Total Observations` = west_n + east_n
  )
})

TableP1_Municipalities <- knitr::kable(
  muni_counts,
  format   = "latex",
  booktabs = TRUE,
  align    = c("c", "c", "c", "c"),
  col.names = c("Bandwidth (km)", "West of the Line", "East of the Line", "Total Observations"),
  caption  = "Number of observations when narrowing the bandwidth"
) %>%
  kableExtra::kable_styling(latex_options = c("hold_position")) %>%
  kableExtra::add_header_above(c("Municipalities" = 4))

print(TableP1_Municipalities)

# --------------------------- #
#         FIGURE P.1.1        #
# --------------------------- #

aggregate_data_12km <- aggregate_data %>%
  dplyr::filter(!is.na(dist_river_combined), dist_river_combined <= 12)

# --------------------------- #
#          PANEL A            #
# --------------------------- #

congress_models_left_12km <- run_models(
  years = congress_years,
  prefix = "spnxleft_congress",
  suffix = "catxleft_congress",
  data = aggregate_data_12km,
  variable_prefix = "left_congress_12km"
)

cp_models_left_12km <- run_models(
  years = cp_years,
  prefix = "spnxleft_cp",
  suffix = "catxleft_cp",
  data = aggregate_data_12km,
  variable_prefix = "left_cp_12km"
)

ols_left_congress_plot_12km <- extract_coefficients(
  model_list = congress_models_left_12km,
  term_of_interest = "east_river",
  model_label = "Left",
  year_labels = congress_years_labels
)

ols_left_cp_plot_12km <- extract_coefficients(
  model_list = cp_models_left_12km,
  term_of_interest = "east_river",
  model_label = "Left",
  year_labels = cp_years_labels
)

all_left_plot_12km <- bind_rows(ols_left_congress_plot_12km, ols_left_cp_plot_12km) %>%
  mutate(term = factor(term, levels = custom_order))

Figure2_PanelA_12km <- ggplot(all_left_plot_12km, aes(estimate, term)) +
  geom_point(size = 2, shape = 19, color = "steelblue", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, size = 0.5, color = "steelblue", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, size = 1.0, color = "steelblue", position = pd) +
  geom_vline(xintercept = 0, linetype = "dotted", color = "gray50", linewidth = 0.3) +
  theme_bw() +
  scale_x_continuous(
    "Estimated effect of living east of the river on voting for the Left",
    breaks = seq(-0.05, 0.20, by = 0.05),
    limits = c(-0.05, 0.16)
  ) +
  scale_y_discrete("") +
  coord_flip() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
  )

print(Figure2_PanelA_12km)

# --------------------------- #
#          PANEL B            #
# --------------------------- #

catalan_congress_models_12km <- run_models(
  years = congress_years,
  prefix = "catxleft_congress",
  suffix = "catxright_congress",
  data = aggregate_data_12km,
  variable_prefix = "catalan_congress_12km"
)

catalan_cp_models_12km <- run_models(
  years = cp_years,
  prefix = "catxleft_cp",
  suffix = "catxright_cp",
  data = aggregate_data_12km,
  variable_prefix = "catalan_cp_12km"
)

ols_cat_cong_12km <- extract_coefficients(
  model_list = catalan_congress_models_12km,
  term_of_interest = "east_river",
  model_label = "Catalan (OLS)",
  year_labels = congress_years_labels
)

ols_cat_cp_12km <- extract_coefficients(
  model_list = catalan_cp_models_12km,
  term_of_interest = "east_river",
  model_label = "Catalan (OLS)",
  year_labels = cp_years_labels
)

cat_all_12km <- bind_rows(ols_cat_cong_12km, ols_cat_cp_12km) %>%
  mutate(term = factor(term, levels = custom_order)) %>%
  arrange(term)

Figure2_PanelB_12km <- ggplot(cat_all_12km, aes(estimate, term)) +
  geom_point(size = 2, shape = 19, color = "#f58231", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, linewidth = 0.5, color = "#f58231", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, linewidth = 1.0, color = "#f58231", position = pd) +
  geom_vline(xintercept = 0, linetype = "dotted", color = "gray50", linewidth = 0.3) +
  scale_x_continuous("Estimated effect of living east of the river on voting for Catalan parties",
                     breaks = seq(-0.05, 0.20, by = 0.05),
                     limits = c(-0.05, 0.16),
                     labels = scales::number_format(accuracy = 0.01)) +
  scale_y_discrete("") +
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 90, vjust = 0.5),
        axis.text.y = element_text(angle = 0, size = 9)) +
  coord_flip()

print(Figure2_PanelB_12km)

# --------------------------- #
#         FIGURE P.1.2        #
# --------------------------- #

aggregate_data_7km <- aggregate_data %>%
  dplyr::filter(!is.na(dist_river_combined), dist_river_combined <= 7)

# --------------------------- #
#          PANEL A            #
# --------------------------- #

congress_models_left_7km <- run_models(
  years = congress_years,
  prefix = "spnxleft_congress",
  suffix = "catxleft_congress",
  data = aggregate_data_7km,
  variable_prefix = "left_congress_7km"
)

cp_models_left_7km <- run_models(
  years = cp_years,
  prefix = "spnxleft_cp",
  suffix = "catxleft_cp",
  data = aggregate_data_7km,
  variable_prefix = "left_cp_7km"
)

ols_left_congress_plot_7km <- extract_coefficients(
  model_list = congress_models_left_7km,
  term_of_interest = "east_river",
  model_label = "Left",
  year_labels = congress_years_labels
)

ols_left_cp_plot_7km <- extract_coefficients(
  model_list = cp_models_left_7km,
  term_of_interest = "east_river",
  model_label = "Left",
  year_labels = cp_years_labels
)

all_left_plot_7km <- bind_rows(ols_left_congress_plot_7km, ols_left_cp_plot_7km) %>%
  mutate(term = factor(term, levels = custom_order))

Figure2_PanelA_7km <- ggplot(all_left_plot_7km, aes(estimate, term)) +
  geom_point(size = 2, shape = 19, color = "steelblue", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, size = 0.5, color = "steelblue", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, size = 1.0, color = "steelblue", position = pd) +
  geom_vline(xintercept = 0, linetype = "dotted", color = "gray50", linewidth = 0.3) +
  theme_bw() +
  scale_x_continuous(
    "Estimated effect of living east of the river on voting for the Left, within 7 km",
    breaks = seq(-0.05, 0.20, by = 0.05),
    limits = c(-0.06, 0.16)
  ) +
  scale_y_discrete("") +
  coord_flip() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
  )

print(Figure2_PanelA_7km)

# --------------------------- #
#          PANEL B            #
# --------------------------- #

catalan_congress_models_7km <- run_models(
  years = congress_years,
  prefix = "catxleft_congress",
  suffix = "catxright_congress",
  data = aggregate_data_7km,
  variable_prefix = "catalan_congress_7km"
)

catalan_cp_models_7km <- run_models(
  years = cp_years,
  prefix = "catxleft_cp",
  suffix = "catxright_cp",
  data = aggregate_data_7km,
  variable_prefix = "catalan_cp_7km"
)

ols_cat_cong_7km <- extract_coefficients(
  model_list = catalan_congress_models_7km,
  term_of_interest = "east_river",
  model_label = "Catalan (OLS)",
  year_labels = congress_years_labels
)

ols_cat_cp_7km <- extract_coefficients(
  model_list = catalan_cp_models_7km,
  term_of_interest = "east_river",
  model_label = "Catalan (OLS)",
  year_labels = cp_years_labels
)

cat_all_7km <- bind_rows(ols_cat_cong_7km, ols_cat_cp_7km) %>%
  mutate(term = factor(term, levels = custom_order)) %>%
  arrange(term)

Figure2_PanelB_7km <- ggplot(cat_all_7km, aes(estimate, term)) +
  geom_point(size = 2, shape = 19, color = "#f58231", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                    xmax = estimate + 1.96 * std.error),
                width = 0, linewidth = 0.5, color = "#f58231", position = pd) +
  geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                    xmax = estimate + 1.645 * std.error),
                width = 0, linewidth = 1.0, color = "#f58231", position = pd) +
  geom_vline(xintercept = 0, linetype = "dotted", color = "gray50", linewidth = 0.3) +
  scale_x_continuous("Estimated effect of living east of the river on voting for Catalan parties, within 7 km",
                     breaks = seq(-0.05, 0.20, by = 0.05),
                     limits = c(-0.07, 0.11),
                     labels = scales::number_format(accuracy = 0.01)) +
  scale_y_discrete("") +
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 90, vjust = 0.5),
        axis.text.y = element_text(angle = 0, size = 9)) +
  coord_flip()

print(Figure2_PanelB_7km)

# --------------------------- #
# --------------------------- #
#         APPENDIX Q          #
# --------------------------- #
# --------------------------- #

# --------------------------- #
#         TABLE Q.2           #
# --------------------------- #

m_pop  <- lm(pop_change ~ east_river, data = aggregate_data)
m_perc <- lm(pop_change_perc ~ east_river, data = aggregate_data)

coef_info <- function(fit, term) {
  s <- summary(fit)$coefficients
  i <- match(term, rownames(s))
  list(b = s[i, 1], se = s[i, 2], p = s[i, 4])
}

stars <- function(p) {
  if (is.na(p)) "" else if (p < 0.01) "**" else if (p < 0.05) "*" else if (p < 0.10) "+" else ""
}

fmt1 <- function(x) formatC(x, format = "f", digits = 1)  # for counts
fmt2 <- function(x) formatC(x, format = "f", digits = 3)  # for percents

b_er_pop   <- coef_info(m_pop,  "east_river")
b_int_pop  <- coef_info(m_pop,  "(Intercept)")
b_er_perc  <- coef_info(m_perc, "east_river")
b_int_perc <- coef_info(m_perc, "(Intercept)")

N1 <- stats::nobs(m_pop)
N2 <- stats::nobs(m_perc)

TableQ2 <- tibble::tibble(
  row = c("East of River", "", "Constant", "", "Observations"),
  `DV: Pop. Change` = c(
    paste0(fmt1(b_er_pop$b),  stars(b_er_pop$p)),
    paste0("(", fmt1(b_er_pop$se), ")"),
    paste0(fmt1(b_int_pop$b), stars(b_int_pop$p)),
    paste0("(", fmt1(b_int_pop$se), ")"),
    N1
  ),
  `DV: Pop. Change (%)` = c(
    paste0(fmt2(b_er_perc$b),  stars(b_er_perc$p)),
    paste0("(", fmt2(b_er_perc$se), ")"),
    paste0(fmt2(b_int_perc$b), stars(b_int_perc$p)),
    paste0("(", fmt2(b_int_perc$se), ")"),
    N2
  )
)

print(TableQ2)
